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ABSTRACT 

We calculate the angular power spectrum of the Cosmic Microwave Background (CMB) temperature 
fluctuations induced by the kinetic Sunyaev-Zcl'dovich (kSZ) effect from the epoch of rcionization 
(EOR). We use detailed A-body+radiative transfer simulations to follow inhomogeneous rcionization 
of the intergalactic medium (IGM). For the first time we take into account the "self-regulation" of 
rcionization: star formation in low-mass dwarf galaxies (10 8 Mq < M < 10 9 Mq) or minihalos 
(10 5 Mq < M < I0 8 Mq) is suppressed if these halos form in the regions that were already ionized 
or Lyman- Werner dissociated. Some previous work suggested that the amplitude of the kSZ power 
spectrum from the EOR can be described by a two-parameter family: the epoch of half ionization and 
the duration of reionization. However, we argue that this picture applies only to simple forms of the 
rcionization history which are roughly symmetric about the half-ionization epoch. In self-rcgulatcd 
rcionization, the universe begins to be ionized early, maintains a low level of ionization for an extended 
period, and then finishes rcionization as soon as high-mass atomically-cooling halos dominate. While 
inclusion of self- regulation affects the amplitude of the kSZ power spectrum only modestly (~ 10%), 
it can change the duration of reionization by a factor of more than two. We conclude that the simple 
two-parameter family does not capture the effect of a physical, yet complex, reionization history caused 
by self-regulation. When added to the post-reionization kSZ contribution, our prediction for the total 
kSZ power spectrum is below the current upper bound from the South Pole Telescope. Therefore, 
the current upper bound on the kSZ effect from the EOR is consistent with our understanding of the 
physics of reionization. 



1. INTRODUCTION 

How was the intergalactic medium (IGM) reionized be- 
fore z = 6? The secondary anisotropy of the cosmic 
microwave background (CMB) at I > 3000 allows us to 
probe the physics of cosmic reion ization via the kinetic 
Sunyaev-Zel'dovich effect (kSZ; iSunvaev k. Zel'dovichl 
119801 ). The temperature of the CMB changes as free 
electrons in ionized gas Compton scatter CMB photons: 
the bulk peculiar velocity of electrons induces Dopplcr 
shifts in the energy of the CMB photons. While the 
spectrum of the CMB remains that of a black body, its 
temperature changes 

Inhomogcncity in the density and velocity of electrons, 
as well as inhomogcncity in ionization fraction, will in- 
duce temperature fluctuations in the CMB, AT/T, given 
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where 7 is the line-of-sight unit vector, v the peculiar 
velocity field, and t the optical depth to Thomson scat- 
tering integrated through the IGM from z = to the 
surface of last scattering at z roc ps 10 3 , where 

dr = c n e (z)<7T( — )dz. (2) 

There are two contributions to the kSZ signal: 

1. Post-reionization contribution. This is the 
contribution from redshifts below z = z ov , where 
z ov is the redshift at which reionization is finished, 
when individual H II bubbles fully overlap with one 
another. While the post-reionization contribution 
depends upon the value of z ov , for which quasar ab- 
sorption spectra suggest z ov ~ 6 — 7, it is not too 
sensitive to the exact value of z ov . We shall not dis- 
cuss this contribution in this paper, but discuss it 
in a subsequent paper (Park et al., in preparation). 

2. Reionization contribution. This is the contri- 
bution from redshifts above z = z ov , where the 
ionization was patchy and incomplete. This con- 
tribution depends not only on z ov , but also on the 
details of the time and spatial variation of inhomo- 
geneous rcionization, which are not yet well con- 
strained; thus, we must explore how predictions 
vary for different models of rcionization. The rcion- 
ization contribution is the main focus of this paper. 
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Modeling the reionization contribution is a challenge, 
as the universe was not ionized homogeneously, but in 
patches. These patches grow over time until they over- 
lap, finishing reionization of the universe. The distribu- 
tion of these patches is determined by non-linear physics: 
non-linear clustering of the sources of ionizing photons; 
non- linear clumping of gas in the IGM; and complex mor- 
phologies of patches resulting from propagation of ion- 
ization fronts in the clumpy IGM. Accurately calculat- 
ing the reionization contribution thus requires numerical 
simulations of cosmological structure formation coupled 
with radiative transfer. 

To model the formation and spatial clustering of the 
sources of ionizing photons, cosmological simulations 
must be performed in a volume large enough to capture 
the crucial spatial variations of this process in a statis- 
tically meaningful way. This requires a volume greater 
than ~ 100 comoving Mpc across, because H II bubbles 
can typically grow as large as ~ 20 comoving Mpc in size. 
These simulations must also have a high enough mass res- 
olution to resolve the formation of the individual galaxies 
which are the sources of ionizing radiation; thus, billions 
of particles are required. The radiative transfer of ioniz- 
ing photons is then calculated on the IGM density and 
velocity fields computed by the cosmological simulation. 

What do current observational data tell us? The South 
Pole Telescope (SPT) experiment has detected an excess 
temperature anisotropy for the CMB on small angular 
scales, which they attribute to the SZ effect. By sub- 
tracting the dominant contribution from the tSZ effect by 
using multiwavelength observations to distinguish it from 
the kSZ effect, the SPT detection yields an upper limit to 
the total kSZ contribution^ The measurements are usu- 
ally reported in terms of the angular power spectrum, 
Ci = TTiErnM 2 - Here, a lm = J d 2 jAT (7)^(7) 
is the coefficient of spherical-harmonics mode, Yz, n , of 
AT. The SPT collaboration reports their measurements 
in terms of the quantity 

D, s i£+M, ,3, 

which we shall compute in this paper. SPT has placed 
an upper bound on the kSZ D t a t / = 3000 of Df% 00 < 
2.8 fxK 2 (jReichardt et all 120121 ). The detection of the 
total SZ effect is complicated by the possible contami- 
nation of the fluctuating signal caused by the cosmic in- 
frared background (CIB) from individual galaxies. The 
kSZ limit loosens to 6.0 fiK 2 when allowance is made 
for a possible correlati on between the thermal Sun yaev- 
Zel'dovich effect (tSZ; IZel'dovich fe Sunvaevl Il969l ) and 
the CIB. Our goal is to see whether these current upper 
bounds are consistent with our models of reionization. 

Following the early analytical calculation done by 
Vishniac for linear density and vel ocity pe r turba- 
tions in a fully ionized medium (jVishniad Il987t 
Uaffe k, Kamionkowskil I1998D . calculations of the kSZ 
effect by cosmic reionization have steadily improved 
over time. Further analytical calculations later incor- 

8 The post-reionization kSZ effect due to the pairwise relative 
motions of galaxy clusters has been d etected by the Atacama Cos- 
mology Telescope (Hand et al. 2012), but its contribution to the 
CMB temperature fluctuation power spectrum has not been de- 
tected yet. 



porated the effects o f inhomogeneous reionization in an 
appro ximate manner (|Gruzinov fe Hu 1998; Sa ntos et all 
12001 . A "semi-numerical" approach was also devel- 
oped by combining the simulated density and velocity 
fields from N-body simulations with a n analytical ansatz 
for tracking the reio nization process (jZahn et al.l 120051 ; 
iMcQuinn et aL|[2005l) . 

Early, pioneering calculations using structure forma- 
tion simulations coupled with radiative t ransfer to model 
inhomogeneous reionization numerically (jGnedin fe Jaffd 
120011 ; ISalvaterra et al.l [2005) underestimated the ampli- 
tude of the kSZ signal, as they used computational boxes 
too small to capture the impact of large-scale velocity 
modes and H II bubbles or an accurate measure of the 
duration of the global EOR. This was demonstrated by 
the first calculations of reionization based on truly large- 
scale (> 100 Mpc) radiative transfer simulations, which 
resolved th e formation of all galact ic halo sources above 
2 x 10 9 M Q (jlliev et al.ll2007HT2008h . These later simula- 
tions demonstrated the importance of a large enough sim- 
ulation volume to capture the effects of long- wavelength 
fluctuations properly. They were also the first to realize 
that it is necessary to correct the kSZ power spectrum 
for the missing velocity power due to the finite box size 
of the simulations. 

For the mass range of galactic halos resolved by these 
simulations, > 10 9 A/q, stars - the sources of reionization 
- were able to form when the primordial composition gas 
inside the halos cooled radiatively by atomic processes 
involving H atoms. They are known as "atomic cool- 
ing halos" to distinguish them from minihalos of mass 
M < 10 8 M Q , with virial temperature T vir < 10 4 K, for 
which star formation is possible only if H2 molecules 
form in sufficient abundance to cool the gas below T v j r 
by rotational-vibrational line excitation. Atomic-cooling 
halos with 10 8 M Q < M < 10 9 M Q also exist and are even 
more abundant than those with M > 10 9 Af Q . These 
low-mass atomic-cooling halos ("LMACHs"), however, 
are prevented from forming stars if they form within 
an ionized patch of the IGM, where the gas pressure of 
the photohcatcd IGM opposes the accretion of baryons 
onto these halos. This "self-regulates" their contribution 
to reionization as the global ionized fraction grows with 
time and more a nd more of these halos are born withi n 
the ionized zones (|Shapiro et alJH99l llliev et al.l[2007al) . 
While the precise value of halo mass which defines the up- 
per edge of this "Jeans-filtered" mass-range is still uncer- 
tain, the high- mass atomic-cooling halos ("HMACHs") 
above ~ 10 9 Af Q are generally free of this suppression. 

To simulate the impact of both LM ACHs and 
HMAC Hs on reionization, it was necessary for llliev et al.l 
(|2007aD to increase their halo mass resolution so as to re- 
solve all the LMACHs, too, by reducing the simulation 
box size to 53 Mpc on a side at fixed N-body particle 
number. This led to the first radiative transfer sim- 
ulations of "self-regulated" reionization, which demon- 
strated the importance of including and then suppress- 
ing the LMA CHs to start reion ization earlier and extend 
its duration (jlliev et al.ll2007aj) . While the end of reion- 
ization is still set by the rapid rise of the HMACHs, in 
that case, when they eventually surpass the saturated 
contribution of the suppressible LMACHs, the effect of 
the LMACHs is to boost the electron-scattering optical 
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depth, r, integrated through the EOR. Such an effect can 
be important for the kSZ fluctuations from the EOR, too, 
but simulating this required us to increase the simulation 
volume again while retaining the high mass resolution re- 
quired to resolve the LMACHs, too. 

Our next generation of simulations involved boxes 163 
Mpc on a side, a volume large enough to predict observ- 
ables like the kSZ effect, but with N-body simulations 
large enough to resolve all halos down to 10 8 M Q and 
incorporate ionization suppression ("Jeans-filtering") of 
the halos of mass between 10 8 M Q and 10 9 M Q CO iev et al.l 
120121) . These smaller -mass halos (LMACHs) are more 
abundant and likely to be more efficient ionizing sources, 
as they may have higher escape fraction and emissivity 
(jlliev et al.l I2012T ). However, as described above, they 
may be suppressed as sources if they form inside ionized 
regions, where ionization heats the gas and makes its 
pressure high enough to resist gravitational collapse into 
such small galaxies. Recently, an additional simulation 
was performed, including this new physics, in an even 
la rger volum e (~ 6 00 Mpc) (Iliev et al. in preparation). 

lAhn et al.1 (|2012f ) expanded the mass range even fur- 
ther by accounting for starlight emitted by minihalos 
(10 5 — 10 8 Mq), as well. In addition to their Jeans-mass 
filtering in ionized regions, they may also be suppressed 
if molecular hydrogen in minihalos is photo-dissociated 
by Lyman- Werner band photons in the UV background 
below 13.6 eV also emitted by the sources of reionization. 
We thus have a simulated model which takes into account 
all the halos down to 10 5 Mq as sources of reionization. 

It is important now to determine if and how the kSZ 
fluctuations from the epoch of reionization arc different 
from the previous predictions when this "self-regulated" 
reionization is taken into account. That is the prime 
focus of th is paper. Som e of ou r results were first sum- 
marized in iShapiro et al.l (120121). 

R ecently. IMesinger et all (12012ft. IZahn et all (|2012h 
and lBattaglia et al.l (|2012aF compared the predicted kSZ 
power spectra from their semi-numerica l calculations, to 
the u pper bounds from the SPT data (Rciehard t et all 
I2012T) . obtaining limits on the epoch and the dura- 
tion of the reionization. Those studies concluded that, 
for a given value of the total Thomson-scattering opti- 
cal depth, the reionization contribution to the kSZ sig- 
nal is mostly sensitive to the duratio n of the reioniza - 
tion denned as Az = Zq 9 % — z?noz, (IZahn et all 120121) 
or z 7 ^ % — (IMesinger et al.l T2012t IBattaglia et all 

l2012al) . IZahn et al.l (|2012f ) claim that the upper bound 



on Dgfooo from the SPT data implies Az < 4 (95% CL) 
for no tSZ-CIB correlation, and Az < 7 (95% CL) for 
the maximum possible tSZ-CIB correlation. However, 
as their methods are based on an analytical ansatz for 
the reionization process, it is necessary to use more self- 
consistent calculations of radiative transfer such as our 
simulation r esults to revisit this issue. We note that 
IZahn et al.l (|2011f ) compared their semi-numerical ap- 
proach to their own numerical simulations using radiative 
transfer, finding an agreement at the level of 50%. 

The remainder of this paper is organized as follows. In 
Section 2, we express the kSZ power spectrum in terms 
of a line-of-sight integral of the transverse momentum 
power spectrum, and show how the transverse momen- 
tum power spectrum is related to the statistics of the 



density and velocity fields of ionized gas. In Section 3, 
we describe the details of the simulations used for our 
study. In Section 4, we present our predictions for the 
kSZ power spectrum and discuss the effects of inhomo- 
geneous reionization as well as of self-regulated reioniza- 
tion. In Section 5, wc compare our results with the recent 
semi-numerical calculations, and show that inclusion of 
self-regulated reionization qualitatively changes the pa- 
rameter dependence of the kSZ power spectrum from 
that without self-regulation. In Section 6, we summarize 
our conclusions. In Appendix A, we give the derivation 
of the kSZ power spectrum written in terms of the trans- 
verse momentum power spectrum. In Appendix B, we 
show how to correct for the missing power due to a finite 
box size of simulations in our method. 

2. BASICS 

2.1. Angular power spectrum of the kSZ effect 

As the Thomson-scattering optical depth, r, is propor- 
tional to the free electron number density, the kSZ effect 
given by Equation ([1]) depends mainly on the specific 
ionized momentum field of the ionized medium, 



xv(l + 8), 



(4) 



henceforth referred to only as "momentum" . Here, \ = 
n e /(nn + 2riHc) is the ionization fraction, and S = (p — 
p)/p is the density contrast of baryons. In general, the 
baryon density is different from the dark matter density, 
especially on scales smaller than the Jeans length. In this 
paper, we shall assume that baryons trace dark matter 
particles, as we are interested in scales bigger than the 
Jeans length of gas at 10 4 K. 
Wc rewrite Equation ([T]) using q as 



AT _ er T n efi 
T ~ c 



e r q-7. 



(5) 



Here, o~t is the Thomson scattering cross section, n e ^ = 
n~H,o + 2nHc,0i is the mean number density of electrons at 
the (fully-ionized) present epoch, and s is the distance 
travelled by photons from a source to the observer in 
comoving units. 

The kSZ angular power spectrum is given b\r| (See Ap- 
pendix [A] for derivation; also see iMa fc Frvl I2002L but 
note that their Equation (4) contains a typo: it is off by 
a factor of (c/H ) 2 ): 



C, = 



ds -2r P q±( k = l / S , S ) ff ~ 



where q_i_(k) = q(k) — fc[q(k) ■ k] is the projection of 
q(k) = J <i 3 x e lk x q(x) on the plane perpendicular to 

the the mode vector k (i.e., q± • k = 0), k = k/|k| 
is a unit vector, and P q± is the power spectrum of q^ 

9 All previous numerical calculations of the kSZ power spectrum 
first created maps using Equation J5]> and then measured C; from 
the two-dimensional Fourier transform of the simulated maps. In 
this paper, we shall use Equation 10 to compute C; using P q± 
measured from three-dimensional simulation boxes at various red- 
shifts, without ever creating maps. While we are the first to apply 
this method to the computation of the kSZ power spectrum, this 
method has been app lied successfully to th e computation of the 
tSZ power spectrum {Rcfrcgicr et al. 2000) as well as to that of 
the power spectrum of anisot ropy of the near infrared background 
IIFernandez et a.l.|[20TOl. [20121) . 
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Fig. 1. — Dimcnsionless power spectra of the curl of the momentum field, k 3 P q± (k)/(2ir 2 ), at z = 9 calculated from the simulation with 
114 h' 1 Mpc in a side. The black solid lines show the raw power spectrum obtained from the TV-body simulation, while the blue lines show 
the power spectrum after being corrected for the missing velocity power due to a finite box size of the simulation. The red lines show the 
missing power added to the black solid lines. The dotted lines show the analytical OV spectrum given in Equation (5J. Left: fully ionized 
case. An excellent agreement between the OV spectrum and the corrected power spectrum shows the validity of our simulation as well as 
that of our method to correct for the missing velocity power. Right: inhomogeneously ionized case, L3. The power spectrum is significantly 
enhanced at k < 1 h Mpc -1 . 



defined by (2n) 3 P q± (k)6 D (k - k') = (q x (k) ■ ql(k')). 
Note that qj_ is often called a transverse (or curl) mode. 
A longitudinal (or gradient) mode is parallel to k and is 
given by q|(k) = fc[q(k) • fc_L ^ 

As we show in Appendix [Al in the small-angle approx- 
imation, the line-of-sight integral cancels out the contri- 
bution from q|| and a half of the power of qx, leaving 
only the remaining half of P q± . This explains a factor of 
two in the denominator of Equation ([6]) . 

Helium atoms are assumed to be singly ionized where 
hydrogen atoms are ionized at least until z ov , the redshift 
which all the H II bubbles overlap to finish the ionization 
of hydrogen atoms due to the similar ionization potential 
of H I and He I. Helium atoms remain singly ionized 
until much later, z ~ 3, after which they are thought 
to be doubly ionized. As we are interested only in the 
epoch of hydrogen-reionization, z > 6, we shall assume 
that the ionized fraction, Xi is given by \ = (0.92)A, 
where X is the hydrogen ionized fraction at each point in 
our radiative transfer simulation: \ saturates at 0.92 in 
fully ionized regions during hydrogen reionization, as 8% 
of the electrons are left bound in singly-ionized helium 
atoms. 

2.2. Power spectrum of the curl of the momentum 

Our goal is to compute the power spectrum of the curl 
of the momentum held, P q± , and evaluate Equation ([6]) 
to obtain C;. 

Assuming that the velocity field stays longitudinal, i.e., 
parallel to k, P q± is given by the second-order term in 

the momentum: q ± = (/ £Mj6(k - k')v(k')) ± . This 



assumption is exact in the linear regime and is approxi- 
mately true in the non-linear regime, as this second-order 
term dominates in the non-linear regime anyway. This 
gives (IMa fc FrvH2"002l) 



P q± (k,z) 



rf 3 V '~> 

' (l-fi' 2 )[P ss (\k-k'\)P vv (k') 



(27T) 



k' 



Ik k 



-P Sv (\k-k'\)P Sv (k') 



(7) 



where fjf = k ■ k' . Here, PssPw term gives a positive con- 
tribution, whereas Ps v Ps v term gives a negative contri- 
bution from the density field correlated with the velocity 
field that does not have a curl component. 

Due to a finite box size of simulations, we must cor- 
rect for the missing velocity power coming from modes 
who se wavelength is lo nger than the size of the simulation 
box (jlliev et al.ll2007bj l We shall describe our correction 
method in Appendix iBl 

At high redshift where the density and velocity fields 
are still in the linear regime, the velocity power spec- 
trum is related to the linear density power spectrum by 
P vv (k) = (af /k) 2 Plf{k), where / = riln(5/rilna and 
a(t) is the Robertson- Walker scale factor. This gives 
the so -called Ostriker-Vishniac (OV) spectrum (jVishniad 
[1981 : 



p°/(M)=a 2 .r 



±^P^(\k-k'\,z)P^(k',z) 

;; fc(fc-2fcV)(l-M' 2 ) 
k' 2 {k 2 + k' -2kk'n')' 



■5 



TABLE 1 

Reionization simulation parameters and global reionization history results 



Label 


9t,h (/t,h) 


9-1, L (/-y,l) 


St.mh (/i,mh£ 




z 10% 


z 90% 




LI 


8.7(10) 


130(150) 




0.080 


13.3 


8.6 


8.3 


L2(XL2) 


1.7(2) 


8.7(10) 




0.058 


9.9 


6.9 


6.7 


L2M1J1 


1.7(2) 


8.7(10) 


5063(1030) 


0.086 


17.4 


6.9 


6.7 


L3 


21.7(25) 






0.070 


10.3 


9.1 


8.4 



a MH efficiencies j 7i mh (/ 7! me) quoted here are for the minimum-mass halo assumed to contribute, 10 Mq, which is roughly comparable 
to the average value for the minihalos integrated over the halo mass function. The efficiency of any MH of a given mass M is obtained 

, 10 5 M , 



simply by multiplying to the quoted 9 7 , M H (/ 7 ,mh) by (- 



-)• 



The OV spectrum provides a useful check of the numeri- 
cal simulation and the way we correct for the missing ve- 
locity. In the left panel of Figure [TJ we show an excellent 
agreement between the OV spectrum and the simulation 
result at z = 9, after correcting for the missing velocity 
power due to a finite box size of the simulation. 

Finally, one can incorporate the effect of inhomogc- 
neous reionization into the equation by replacing 5 in 
Equation (7} by x(l + 5): 



P q± (k,z) 



afk 

W) 



3 (1 - M ' 2 ) [P x(1+5) , x(1+5) (|k - k'\)P vv (k') 



|k-k 



7r^(i+5),,(|k-k'|)P x(1+5)it ,(fc') 



(9) 



Note that we do not use this equation to compute P q± , 
but compute P q± directly from the simulation. How- 
ever, we use this equation to estimate and correct for the 
missing power due to a finite box size of the simulation 
as described in Appendix IBl We then use the corrected 
P q± in Equation ([6]) to compute the angular power spec- 
trum. As shown in the right panel of Figure [U the effect 
of reionization inhomogeneity substantially boosts the 
power spectrum relative to the homogeneously-ionized 
case, while correcting for the missing velocity power of 
the finite simulation volume boosts it even further. 

3. REIONIZATION SIMULATION 
3.1. Basic simulation parameters 

The simulations that we shall use in this paper con- 
sist of two parts: (1) cosmological iV-body simulations 
of coll isionless particles using the "CubeP 3 M" A^-body 
code (jHarnois-Deraps et al.l 120121) ; and (2) radiative- 
transfer of H-ionizing photons in the density and source 
fields created from this N-body simulation results using 
the "C 2 -Ray" (Conse rvative, Causal Ray-tracing) code 
(jMellema et al.l 120061 ). The details of the simulations 
that w e sh all use in thi s pape r are described in llliev et al.l 
(12011) and lAhn et all (12011 . 

Unless specified otherwise, the reionization simulations 
are run on the density and source fields from the same 
N-body results with 3072 3 particles in a comoving box 
of 114 /i -1 Mpc on a side. Halos are identified down 
to 10 s Mq with at least 20 particles, using a spherical 
overdensity halo finder with overdensity of 178 times the 
mean cosmic density. One of the models uses another N- 
body simulation with a larger box of 425 Mpc, with 
5488 3 particles, resolving halos down to 10 9 Mq. The 
background cosmology is based on the WMAP 5-year 



data combined with constraints from baryonic acoustic 
oscillations and high-redshift Type la supernovae (J7m = 
0.27, » A = 0.73 h = 0.7, Q b = 0.044, <j 8 = 0.8, n s = 0.96; 
iKomatsu eraTll2009[ ). 

For the 114 /i" 1 Mpc, we then calculate the IGM 
density field from the particle data with halos excluded 
adaptively-smoothed on to a 256 3 radiative-transfer grid 
in order to generate ionization maps using the C 2 -Ray 
code. Therefore, the final physical length resolution of 
the reionization models is d co n = 0.45 hr 1 Mpc. The 
highest Z-modc that we can calculate from the simulation 
is given by iumit = fcNy q s(-z ov ), where fc Nyq = 7r/(2d ceU ) 
is the Nyquist frequency, and s(z ov ) is the comoving 
distance out to the end of reionization. For example, 
z ov = 6.6 gives /n mit = 22000. 

The new simulations also incorporate the effects of 
even smaller halos in 10 5 Mq < M < 10 s Mq, using a 
sub-grid prescription calibrated by smaller-box N-body 
simulations with higher-re solution having 1 728 3 particles 
in a box of 6.3 h' 1 Mpc (|Ahn et al.ll2012D . Specifically, 
we find that there is a correlation between the num- 
ber of these small-mass halos in each cell and the total 
matter density averaged over that cell, with cells of size 
0.45 hr 1 Mpc, which coincides with the size of the radia- 
tive transfer cells in our 114 hr 1 Mpc C 2 -ray simulations. 
We then use this correlation to calculate the number of 
small-mass halos in each of the radiative-transfer cells in 
our 114 hr 1 Mpc simulations. 

For our most recent simulation, in a box 425 Mpc 
on a side, the RT grid has 504 3 cells, so d cc \\ = 
0.84 hr 1 Mpc, slightly larger than that for the other 
simulations, and /n m it ~ 12000. In this larger-box sim- 
ulation, low-mass halos between 10 8 and 10 9 Mq are 
included by a subgrid model like that described above 
for MHs. 



3.2. Varying physics of reionization 

What kind of sources are responsible for reionization? 
In this section, wc consider a set of reionization simula- 
tions based on source models of increasing sophistication 
from the one with only high-mass sources to the one with 
all kinds of sources down to least massive halos in our 
models. 

For each halo identified in our simulation, we calcu- 
late the number of ionizing photons which escape from 
it into the IGM per unit time, N~, which is assumed to 
be proportional to the halo mass, M: 



N~ 



f 7 MQ b 
At Clomp 



(10) 



6 



where m p is the proton mass, At is the duration of each 
star-forming episode (i.e. which corresponds in prac- 
tice to the radiative transfer simulation time-step), and 
/ 7 = /esc/*-^V* is the number of ionizing photons pro- 
duced and released by the halo over the lifetime of the 
stars which form inside it in this time step, per halo atom, 
if /* is the fraction of the halo atoms which form stars 
during this burst, / esc is the fraction of the ionizing pho- 
tons produced by these stars which escapes into the IGM 
and the integrated number of ionizing photons released 
over their lifetime per stellar atom is given by AT*. The 
latter parameter depends on the assumed IMF for the 
stellar population and can range from ~ 4, 000 (e.g. for 
Pop II stars with a Salpctcr IMF) to - 100, 000 (e.g. for 
a top-heavy IMF of Pop III stars). Halos were assigned 
different efficiencies according to their mass, grouped ac- 
cording to whether their mass was above ("HMACHs") 
or below ("LMACHs") 10 9 M© (but above 10 8 M Q , the 
minimum resolved halo mass). Low-mass sources are as- 
sumed to be suppressed within ionized regions (for ion- 
ization fraction higher than 10%), throug h Jeans-mass 
filtering, as discussed in llliev et al.l (|2007aD . 

In addition to the source efficiency parameter, / 7 , we 
also define a slightly different factor, g 7 , that is given by 

* - A (^) <") 

where At is the time between two snapshots from the 
N-body simulation. The new factor c/ 7 reflects the fact 
that a given halo has a luminosity which depends on 
the ratio of / 7 to At, so <? 7 has the advantage that it is 
independent of the length of the time interval between 
the density slices, and as such it allows a direct com- 
parison between runs with different At. For the reader's 
convenience, we listed the values of both parameters in 
Table [TJ The specific numerical values of the efficiency 
parameters are strongly dependent on the background 
cosmology adopted and the minimum source halo mass. 
Therefore, parameter values for simulations based on dif- 
ferent underlying cosmology and halo mass resolution 
should not be compared directly, but require cosmol- 
ogy and resolution-dependent conversion coefficients to 
achieve the same reionization history. 

3.2.1. HMACHs-only model 

In our simplest model (labeled as L3; see TableQ]for the 
details. Note that "L" stands for a "large volume"), we 
only use HMACHs as the sources of reionization. These 
sources are defined as the halos with M > 2.2 x 10 9 M Q 
for L3; and with M > 10 9 M Q for the other configura- 
tions. These sources are believed to form stars even when 
immersed in ionized regions, due to the fact that their 
gravitational potential wells are deep enough to overcome 
Jeans-mass filtering. 

3.2.2. HMACHs+LMACHs models 

What about smaller-mass halos? LMACHs are more 
abundant; however, if they form inside the regions that 
have already been ionized, they would not act as sources 
of ionizing photons. This is because ionization heats the 
gas and makes its press ure too high for th e gas to collapse 
into such small halos (|Iliev et al.ll2007al and references 
therein) . 



When we include LAMCHs and account for this "self- 
regulation" of reionization, we give LMACHs a higher 
efficiency, g 7 , than for HMACHs, as presumably it is 
easier for ionizing photons to escape from LMACHs than 
from HMACHs, and Pop III stars with a top-heavy ini- 
tial mass function (IMF) , which are capable of producing 
more ionizing photons than Pop II stars with a Salpeter 
IMF, are more likely to form in LMACHs. 

There are two cases which have both HMACHs and 
LMACHs, and we shall call them LI and L2. For LI, the 
efficiency parameter, c/ 7 , is chosen such that the overlap 
redshift, z ov = 8.3, is similar to that of L3, z ov = 8.4 
(see Table [Q. For L2, <? 7 is chosen such that z ov is be- 
tween 6 and 7, as suggested by the quasar absorption line 
observations. 

For L2, we have another run with a much larger 
volume (425 h~ x Mpc) with 504 3 of radiative-transfer 
grids. Although it does not resolve LMACHs, we include 
LMACHs as a sub-grid model using correlation between 
average density of radiative transfer cells and number 
density of LMACHs similarly to how lAhn et al.l (|2012f ) 
included MHs in the simulation (Iliev et al. and Ahn et 
al. in preparation). This run gives Ziimit ~ 12000. We 
shall call this configuration "XL2" , as the volume for this 
run is bigger (hence the name, XL) than those runs with 
"L." This run will be used to check our method to correct 
for the missing velocity power. 

3.2.3. HMACHs+LMACHs+MHs model 

What about even smaller-mass sources? Gas in halos 
of masses between 10 5 M@ and 10 s M is thought to 
cool via rotational and vibrational transitions of hydro- 
gen molecules and form stars, until hydrogen molecules 
are dissociated by Lyman- Werner photons in the UV 
background from other sources (see lAhn et al.ll2"012t and 
references therein). 

The MHs form earlier than LMACHs or HMACHs, 
and thus can start reionization of the universe earlier. 
However, as the star formation in MHs is vulnerable to 
Lyman- Werner photons, it gets suppressed wherever the 
intensity of the LW background rises above the thresh- 
old for suppression, locally at first, and eventually glob- 
ally. This adds another kind of "self-regulation" to the 
reionization history, with an even more extended phase 
of low-level ionizati on before MHs are eventually sup- 
pressed completely (jAhn et a l. 2012). 

The effects of MHs have been added to L2 bv lAhn et al.l 
HoH), and we take one of the cases simulated there. 
L2M1J1, as our fiducial case with MHs. See Table Q] for 
the efficiency of MHs. "M" denotes the mass spectrum 
of Pop III stars in MHs, and "J" the threshold intensity 
of the Lyman- Werner photon background, above which 
the star formation in MHs is suppressed. In L2M1J1, 
each halo is assumed to host one Pop III star with mass 
of 300 Mq, and the assumed LW threshold is JLw,th = 
0.1 x 10~ 21 ergs _1 cm _2 sr _1 . 

4. RESULTS 

Before presenting and discussing our predictions for 
the kSZ power spectrum, let us briefly comment on the 
global ionization history of the universe, which is the key 
to understanding the difference between our results and 
the previous ones. For more detailed discussion on the 
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Fig. 2. — Cuts through the 7V-body+Radiative Transfer simulations used in this work. See Table [Tl for the parameters of models LI, L2, 
L2M1J1, and L3. While these runs have the box size of 114 h~ Mpc, the model XL2 has the box size of 425 /i -1 Mpc and has the same 
model parameters as the model L2. Each panel shows the matter density distribution multiplied by spatially-varying ionization fractions. 
For example, it just shows the matter density when a given region is fully ionized, while it shows nothing (i.e., white) when a given region 
is fully neutral. The density fields are color-coded such that overdense regions are red and undcrdense regions are blue. We create this 
figure by interpolating between adjacent snapshots at a given lookback time. The length scale is linear in the co-moving units. The x-axis 
shows redshifts, while the jy-axis shows h~ 1 Mpc. 



effects of self-regu lation, see llliev et al.l (|2007al [2012t) and 
lAhn et all IpOia ). 

Figure [5] shows how the reionization proceeds in 
our simulation boxes, while Figure [3] shows the mass- 
averaged ionization fraction of the universe as a function 
of redshift. Both figures show that inclusion of low- mass 
halos (LMACHs and MHs), which are self- regulated, sig- 
nificantly extends the ionization history of the universe 
toward higher redshift. Let us compare LI and L3. As 
LMACHs form earlier, the universe begins to be ionized 
earlier in LI than in L3. However, the universe does not 
get reionized quickly but keeps a low level of ionization 
for an extended period due to self-regulation of sources. 
Only after HMACHs start to dominate, at z ~ 10, does 
reionization proceed rapidly and finishes soon thereafter. 
In L3, with no LMACHs, by contrast, reionization pro- 
ceeds rapidly from beginning to end because the abun- 
dance of HMACHs, the only sources, grows exponentially 
without any suppression effects to sclf-regulate them. 
When MHs arc included (L2M1J2), the universe begins 
to be ionized even earlier than the cases with HMACHs 
and LMACHs, and keeps a low- level ionization for a 



longer period. 

These physically motivated yet somewhat complex 
reionization histories were not considered in any of the 
previous calculations of the kSZ power spectrum. In 
this section, we show that it is these new features in the 
reionization history that invalidate simple two-parameter 
descriptions of the amplitude of the kSZ power spec- 
trum proposed by the previous study (iZahn et al.1 120121 ; 



iMesinger et al.ll2012l ; iBattaglia et alil2012al ) 

4.1. Impact of Inhomogeneous Reionization 

First, it is useful to understand how important it is 
to include inhomogeneity (or patchiness) of reionization 
when computing the kSZ power spectrum. In order to 
see this, we create a homogeneous version of L3 ("L3- 
homogeneous" ) , in which we wipe out inhomogeneity of 
reionization by replacing the ionization fraction, Xi with 
its global average, x ( see Figure [3]). This then gives the 
transverse momentum power spectrum as P q± = x 2 P^ , 
where P^ v is the OV spectrum given by Equation ([8]). 
We remind reader that, on the scales of interest to us in 
this power spectrum (k < 1 h Mpc -1 ), the degree of non- 
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Fig. 3. — The global mean ionization history of our models (see Table [T] for the parameters of models). The mass-averaged hydrogen 
ionization fraction, X, is plotted against z. Note how self-regulation results in an extended period of low-level ionization by com paring the 
case without self-regulation (L3 = HMACHs only) and that with self-regulation (LI = HMACHs + LMACHs) l|Iliev et al.H2012l ). A further 
exten sion occurs when MH sources are included, as well (i.e. compare L2 = HMACHs + LMACHs and L2M1J1 = L2 + MHs) l|Ahn et al.l 
12013) . 



linearity of the underlying density and velocity fields of 
the IGM is small enough that we can well approximate 
the kSZ power spectrum for this "homogeneous" ioniza- 
tion case by the assumption of linear perturbations in- 
herent in Equation ([5]) (see Section |2"T21 and the left panel 
of Figure [T]). We use this momentum power spectrum in 
Equation ([6]) to obtain the kSZ power spectrum for "L3- 
homogencous." Thus, "L3" and "L3-homogeneous" have 
exactly the same average reionization history, while spa- 
tial fluctuations of ionization fraction arc included only in 
L3. We find that L3 yields an order-of- magnitude larger 
power spectrum than L3-homogene ous that is consistent 
with findings in llliev et all (|2007bD (see Figure |4]). 

In order to see the effect of inhomogeneous reioniza- 
tion on the kSZ power spectrum in more detail, we show 
the contribution from a given comoving distance to the 
kSZ power spectrum at / = 3000, dC§| 000 /ds, in Fig- 
ure [5J While both L3 and L3-homogeneous converge to 
the same dCf^ 000 /ds after the universe becomes fully 
ionized, we find a clear enhancement of the power when 
the ionization fraction is less than unity, z > z ov = 8.4. 
The maximum contribution occurs when the universe is 
half ionized. One can see this visually in the middle (L3) 
and bottom (L3-homogeneous) panels of Figure [5j L3 is 
clearly more patchy than L3-homogeneous. 

The angular scale for I = 3000 roughly corresponds to 
the co-moving length of 15 h~ x Mpc during the reioniza- 
tion era (z ~ 10). The contribution to the kSZ power 
spectrum continues to grow until the typical comoving 
size of ionized bubbles reaches 15 Mpc. In our mod- 



els, this occurs when the universe is half ionized. After 
this epoch bubbles grow bigger than 15 h^ 1 Mpc, and 
thus the ionization field is no longer patchy on the scale 
of 15 h^ 1 Mpc. This explains why the contribution to 
the kSZ power spectrum at / = 3000 decreases after the 
half-ionization epoch. (By the same token, a plot like 
that for the inhomogeneous case L3 in Figure [5] but for 
/ > 3000 would look similar but with the peak shifted to 
higher z, when ionized patches were smaller-scale.) 

4.2. Impact of LMACHs 

How does the presence of LMACHs and self-regulation 
affect the kSZ power spectrum? To answer this we com- 
pare LI and L3, which arc mostly similar except that LI 
has low-mass halos (10 s M < M < 2.2 x 10 9 M©) with 
most of them being LMACHs. While they finish reion- 
ization at nearly the same redshift, LI begins ionization 
earlier due to LMACHs and gives an extended period of 
low ionization due to self- regulation (see Figure [3]). 

Figure S] shows that LI and L3 give similar kSZ power 
spectra at I < 3000, while at higher multipoles LI be- 
comes significantly greater than L3. This is because 
there are numerous ionized bubbles created by LMACHs 
at high redshifts, which give significant contributions to 
the small-scale kSZ power spectrum. Although it would 
be a challenge for current surveys, future measurements 
of Df sz with 10% accuracy over a wide range of multi- 
poles can distinguish between the predictions of LI and 
L3, shedding light on the roles of LMACHs during the 
reionization. 
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Fig. 4. — Predicted kSZ power spectra, Df° z , for the models discussed in this work (see Tablc[T]for the parameters of models). The box 
size of LI, L2, L2M1J1 and L3 is 114 h _1 Mpc, while that of XL2 is 425 h -1 Mpc. The model parameters of XL2 are the same as those 
of L2, and thus XL2 provides a useful check of the way we correct for the missing velocity power in 114 h _1 Mpc-box simulations (see 
Appendix B for details). The primary CMB power spectrum is also shown. 
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Fig. 5. — Left: The top panel shows the contribution from a given comoving distance to the kSZ power spectrum at I = 3000, dCY^Q 00 /ds. 
The solid line with a peak shows L3, the dashed line shows L3- homogeneous, and the nearly-horizontal solid line shows the fully-ionized 
case. The middle panel is the same as the bottom panel of Figure[2] The bottom panel shows L 3- homogeneous, i.e., the density distribution 
multiplied by the average ionization fraction. Right: A snapshot of L3 at z = 9.3, which gives the maximum contribution to the kSZ power 
spectrum at I = 3000. 
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Fig. 6. — Same as the left panel of Figure[5] but for comparing LI (bottom panel) and L3 (middle panel). See Table[T]for the parameters 
of LI and L3. 
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Fig. 7. — Cumulative reionization kSZ power spectrum at I = 3000 as a function of the maximum redshift (Left) and the mean ionization 
fraction (Right). 
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Fig. 8. — Same as the left panel of Figure [5] but for comparing L2 (bottom panel) and L2M1J1 (middle panel). See Table [T] for the 
parameters of L2 and L2M1J1. 



We compare the contributions from a given comov- 
ing distance to the kSZ power spectrum at I = 3000, 
dCf=3 000 /ds, for LI and L3 in Figure^ As expected, LI 
has larger contributions at higher rcdshifts (z > 10) due 
to LMACHs. On the other hand, L3 has larger contribu- 
tions at lower redshifts (z < 10), as it is more patchy due 
to the absence of smaller bubbles around LMACHs (see 
the middle (L3) and bottom (LI) panels of FigureE]). In 
LI, bubbles around LMACHs do not grow much because 
of self-regulation. 

In the left panel of Figure [7j we show the cumulative 
contributions to the kSZ power spectrum at I — 3000 be- 
low a given maximum redshift, z. This also shows that LI 
receives larger contributions from higher redshifts than 
L3: 20% of the total power in LI comes from z > 11, 
while only a few percent of the total power in L3 comes 
from z > 11. Similarly, the right panel of Figure [7J shows 
that 20% of the total power in LI comes from when the 
ionization fraction is less than 0.25, which is consistent 
with the ionization history above z = 11 shown in Fig- 
ure [3] This extended tail has important implications for 
the interpretation of the kSZ power spectrum, as we shall 
discuss in Section [5j 

4.3. Impact of Mini Halos 

What about MHs? We compare L2 and L2M1J1, 
which have the same efficiency parameters for HMACHs 
and LMACHs, but only L2M1J1 considering MHs. While 
L2 and L2M1J1 finish reionization at almost the same 
redshift, L2M1J1 begins ionization much earlier due to 
MHs and gives a significantly more extended period of 
low ionization due to self- regulation (see Figure [3]). 

Figure 0] shows that L2 and L2M1 Jl give similar kSZ 
power spectra at I < 5000, while at higher multipolcs 
L2M1J1 becomes greater than L2. The reason is the 
same as that for LI versus L3: there are numerous ion- 
ized bubbles created by MHs at high redshifts, which 
contribute to the small-scale kSZ power spectrum. 

While L2M1J1 begins reionization much earlier and 



thus has more contribution from high redshifts to the 
kSZ power spectrum, the actual magnitude of the high- 
rcdshift contribution is modest. This is because of self- 
regulation: self-regulation prevents bubbles around MHs 
from growing, and thus we end up having numerous small 
bubbles filling space nearly uniformly. This results in a 
lesser degree of patchincss, hence a modest contribution 
to the kSZ power spectrum at I = 3000. One can see this 
visually in the middle (L2M1J1) and bottom (L2) panels 
of Figure [5J As a result, the situation is similar to that 
for LI versus L3: 20% of the total power at I = 3000 in 
L2M1J1 comes from z > 10, while only 5% of the total 
power in L2 comes from z > 10. 

It is interesting that all the models with self-regulation 
(LI, L2, and L2M1J1) lie on top of each other when the 
cumulative contribution is shown as a function of the 
mean ionization fraction (see the right panel of Figure[7j, 
whereas the model that docs not have self-regulation (L3) 
is a clear outlier. Whether this is merely a coincidence or 
a unique feature of self-regulation is unclear due to the 
limited number of samples. 

5. SPOT CHECKING THE PREVIOUS CONSTRAINTS ON 
THE DURATION OF REIONIZATION: MORE 
EXTENDED HISTORIES CAN GIVE SIMILAR KSZ 
SIGNALS 

What determines the amplitude of the kSZ power 
spectrum? Rec ent studies using semi-numerical reion- 
ization models (IZahn et all 120121: IMesinger et~aTI 120121: 
IBattaglia et~alll2012al) claim that the amplitude of the 
kSZ power spectrum at I = 3000 can be described by 
a two-parameter family: the redshift of half-ionization, 
z 50%i an d the duration of reion ization defined as Az = 
fQ 3 % - z? M (IZahn et all I2012D or A z = z 7 5% - z 25% 
(jMesinger et all 120121 : IBattaglia et afll2012al) . None of 
these studies included the effects of self-regulated reion- 
ization, and thus the reionization histories explored in 
these studies are roughly symmetric about the epoch of 
half-ionization. 

Figure 2 of IZahn et all (|2012f) shows that the kSZ 



12 



TABLE 2 

Global reionization history and kSZ signal 



Label 


Z 5Q% 


z 99% - z 20% 


z 75% ~ z 25% 


^JSnon^ > 5 - 5 ) 


D fS<mo(< z ov)l 




r-jkbZi. total 


LI 


9.5 


3.2 


2.2 


1.27 


1.94 


0.83 


2.77 


L2 


7.6 


2.1 


1.4 


0.87 


1.69 


0.66 


2.35 


L2M1J1 


7.7 


6.5 


2.1 


0.90 


1.69 


0.69 


2.38 


L3 


9.1 


1.3 


0.9 


1.20 


1.96 


0.75 


2.71 



a From the scaling relation of Shaw ct al. (2012). 

power spectrum at / = 3000 increases by a factor of two 
as the duration o f reionization increase s from Az = 2 
to 4. Figure 10 of iMesinger et all (|2012[) shows that, for 
a half-ionization redshift of 250% = 9, the kSZ power 
spectrum at I = 3000 increases by a factor of 1.4 as the 
duration of reionization increases from Az = 1.3 to 2.6. 
The former gives a scaling of Df^ 000 oc (z 9g% - z 20% ), 

whereas the latter gives Df^ 000 oc (275% — z 25 %) - 5 , for 
a fixed half-ionization redshift. More recently, using a 
new semi-numerical method based on a correlation be- 
tween the smoothed density field and the redshift-of- 
re ionization field found from radiation- hydro simula tions 
of iBattaglia et all (|2012b[ ). IBattaglia et all (|2012afl cal- 
culate the kSZ power spectrum coming from z > 5.5 and 
obtain the following scaling relation: 



D^> 55 = 2.02 /zK 2 



1 + z 
11 



- 0.12 



A: 



0.47 



1.05, 

(12) 

where Az = z 75 % — z 25 % and z is the mean value of 
the redshift-of-reionization field, which is approximately 
equal to the half-ionization redshift, 250%. 

Our predictions for Df^ QQ0 are summarized in Tabled 
Among the models we have explored in this paper, L3 
(which contains only HMACHs and does not have self- 
regulation) closely matches the scenarios explored in the 
above studies. Using z 50 % = 9.1 and z 75 o/ — z 25 % = 0.9 
we find for L3, Equation jT2j gives D^ z > 5 5 = 1.5 /xK 2 . 
This is in a reasonable agreement with our result PI 

£,kSZ,z>5.5 , „ K 2 

However, the above formula significantly overestimates 
the amplitude of the kSZ power spectrum for LI: Equa- 
tion (fT2"j) gives -D; = 3ooo 5 ' 5 = 2.4 ^K 2 , whereas we find 

^kSZ,z>5.5 _ 1 q , llr 2 



(=3000 



1.3 ptK . In other words, despite the fact 
that Ll has a significantly more extended duration of 
reionization than L3 (by a factor of more than two), 
z 75% ~ z 25% = 2.2, the amplitude of the kSZ power spec- 
trum increases only by 8%. Similarly, Equation (|12|) gi 



kSZ,z>5.5 
(=3000 



1.5 and 1.9 fiK 2 for L2 and L2M1J1, respec- 
tively, whereas we find 0.9 fiK 2 for both cases. Therefore, 
we conclude that Equation (fl~2"j) is valid only for simple 
scenarios where the reionization history is roughly sym- 
metric about the half-ionization redshift, but is invalid 
when self-regulation is in clud ed. Similar conc l usion s ap- 
ply to lZahn et all (|2012D and IMesinger et al.l (|2012t) . 
Our results show that self-regulation makes the dura- 



10 In order to compute £^3000 8 S 
from 2 between 5.5 and z, 
P° v , and add it to D ; k |f 00() (> 



of Table [2] 



we calculate the contribution 
le fully-ionized formula, P q± = 
) shown in the seventh column 



tion of reionization significantly more extended without 
changing the amplitude of the kSZ power spectrum very 
much. In other words, an extended period of low-level 
ionization in z > z 50 % does not make much contribution 
to the kSZ power spectrum at I = 3000. 

6. CONCLUSION 

In this paper, using the state-of-the-art reionization 
simulations incorporating the effects of self -regulated 
reionization (|Iliev et~aT1l20TirAhn et al.N20H . we have 
computed the power spectrum of the kSZ effect from the 
EOR. Unlike the previous work which created maps and 
computed two-dimensional Fourier transforms from the 
maps, we have computed the kSZ power spectrum from a 
line-of-sight integral of the transverse momentum power 
spectrum of ionized gas. We present a method to sta- 
tistically correct for the missing velocity power in Ap- 
pendix^ and verify the accuracy of our method by com- 
paring the results from large- (425 Mpc / h) and small-box 
(114 Mpc/h) simulations. 

We find that the kSZ power spectrum is a sensitive 
probe of patchiness of reionization: patchiness increases 
the amplitude of the kSZ power spectrum by an order 
of magnitude. The maximum contribution occurs when 
the angular sizes of ionized bubbles arc close to those 
corresponding to a given multipole. 

While inclusion of small-mass halos such as LMACHs 
and MHs makes the beginning of reionization earlier, 
self-regulatio n significantly slows down the progress o f 
reionization (|Iliev et al.l l2007aL l20ll lAhn et alJ l20Tl . 
This results in an extended period of low-level ioniza- 
tion before more massive HMACHs dominate and finish 
reionization. We find that such an extended period of 
low-level ionization does not, however, make much of a 
contribution to the kSZ power spectrum at / = 3000: 
-C^Jfooo crian g es only by ~ 10% despite the fact that 
the duration of reionization increases by a factor of more 
than two. 

Our results qualitatively change the conclusions 
reached by the previous work w hich did not in- 
clude self-regul ation. Recent work flZahn et aLl 120121 : 
IMesinger etaLI 120121: IBattaglia et al] I2012al) ~assumes 
that -D; k Jf oo can he adequately parameterized by the 
redshift of half-ionization, 2:50%, and the duration of 
reionization, Az. While our result for the simplest model 
of reionization without self-regulation (L3) agr ees with 
the scaling formula of IBattaglia et al.l (|2012af) (Equa- 
tion IT2"j) , our results for the models with self- regulation 
do not agree with it: specifically, the amplitude of the 
kSZ effect is no longer correlated well with the dura- 
tion of the reionization. This is because self-regulation 
gives an extended period of low-level reionization only 
for z > z 2 q%, while the simple models adopted by these 
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other treatments have a roughly symmetric reionization 
history about z = z 50 % , for which a longer duration thus 
implies a longer period of patchy state with a signifi- 
cant ionization across z = z 50 %. Therefore, a more accu- 
rate scaling formula is required to take into account the 
asymmetric reionization history typical of self-regulated 
reionization. 

Going beyond I = 3000, we find that LMACHs 
and MHs do have a considerable impact on the kSZ 
power spectrum on smaller angular scales. For example, 
Dfjfoooo is boosted by 60% and 25% when LMACHs and 
MHs arc included, respectively. Even though measure- 
ments of the kSZ power spectrum at I > 3000 would be a 
challenge for the moment due to contamination by extra- 
galactic point sources and tSZ, future multi-wavelength 
observations may allow us to determine the kSZ power 
spectrum from the EOR over a wide range of multipoles. 
Such measurements will provide us with valuable addi- 
tional information on the nature of the ionizing sources 
and the history of reionization. 

How do our calculations compare with these current 
observational constraints? In order to obtain the total 
kSZ signal from both reionization and post-reionization 
contributions, we take the "CSF" (c ooling and s t ar for - 
mation) post-reionization model of iShaw 
that approximately incorporates the Jeans-filtering of 
P q± due to shock heating in halos and in the IGM. The 
post-reionization kSZ signal computed from their scal- 



ing relation and the total kSZ signal (i.e., the sum our 
reionization calculation and their post-reionization cal- 
culation) are shown in the sixth and seventh columns of 
Table [3J respectively. We find that all of our predictions 
are consistent with the 95% CL upper bound on the total 
signa l from SPT, D^°* al < 2.8 fiK 2 (|Reichardt et al.1 
[2012T ). Therefore, we conclude that the current data are 
consistent with our understanding of the physics of reion- 
ization. 
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APPENDIX 

DERIVATION OF THE POWER SPECTRUM OF THE KSZ EFFECT 

Suppression of longitudinal modes 

An important observation of the nature of kSZ is that it is given by the transverse (vector-mode or spin-1) momentum 
field, and the longitudinal contribution is suppressed. To show this, we Fourier transform Equation ((5|): 



AT _ o- T n e fi 
T c 



ds 



d 3 k 
(2^ 



[7- q(k, s)]e 



— zk- (57) 



(Al) 



Decomposing the momentum vector in Fourier space, q, into the longitudinal component, q\\ = q-fc, and the transverse 
component, gj_ = |q — fc(q ■ fc)|, we obtain 

d 3 k 



AT _ o- T n e fl 
T ~ c 



ds 



(2tt)s 



xq\\ (k, , 



cos 



<^)(l-x 2 ) 1/2 9J-(k,s) 



(A2) 



where x = k ■ 7, and <j)q and are the angles that k makes with q and 7, respectively. 

If the factor e tksx oscillates much more rapidly than the other quantities, the integral over s will be small due to 
cancellation. Recalling that a(s), t(s), and q all vary over the Hubble length scale, kx should be much smaller than 
H/c in order to avoid the cancellation. Namely, either the wavelength should be longer than the Hubble length, i.e., 
k < H/c, or the mode should be nearly perpendicular to the line-of-sight direction, i.e., x ss 0. The former docs 
not contribute much because the amplitude of such a long-wavelength mode is small. Thus, only the modes that are 
perpendicular to the line-of-sight direction, x ~ 0, have a chance to contribute to the kSZ signal. 

However, in this configuration, the longitudinal component of the momentum field is also perpendicular to the line- 
of-sight, and vanishes when taken a dot-product with the line-of-sight, i.e., xq\\ ~ 0. Therefore, only the transverse 
mode survives in the integral, giving 



AT 

— (7) 



ds 



d 3 k 
(2^)3 



<^)(l-x 2 ) 1/2 9±(k,s)< 



— iks 



(A3) 



Angular Power Spectrum 

Here, we follow steps similar to those in Chapter 7.3 of IWcinberd (|2008f ) to derive the angular power spectrum of 
CMB fluctuations induced by the kSZ effect. 
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Spherical harmonic decomposition of Equation (|A3[) gives 

/AT 
d 2 lYril) — (7) 



where 



I d ^y rW / ^e-y ^cos(^-^)(l-x 2 )^^( k , s ) e -^ 
cr T n ei0 /" d 3 k 



jd^Yril) J ^ I e"^cos(^-^)(l~ a; 2 ) 1 / 2 ^(k, S )e-^ 

d 2 7^""(7) / -^ I e^cos(^-^)(l- a; 2 ) 1 / 2 ^(k, S ) 
x 4tt ^(-i) L iL (MF L M (7)^l M * (*) ■ ( A5) 

We first choose a convenient coordinate system in which the z-direction lies on that of the mode vector, i.e., k = z, 
and the azimuthal direction is the same as the direction of the momentum vector, i.e., <fiq = 0. In this case, Y^f*{k) 

simplifies to Y^*{z) = S M0 J^^, giving 



T~ / ^2 e ^(k, S )^(- l ) L V2TTTjL(M J d 2 7 Ylil) [Yr\l) - Y?(V] Yr*(% (A6) 



where 9 and </> = 4"y determine the line-of-sight vector as 7 = (cos 9 sin </>, sin 9 sin </>, cos ( 
The integral over 7 can be computed using, 



d 2 l Yt\i) 1T( 7 ) Yr*{l) = y ^^l+r) 1 ^ ^^' M; m ' ~ Ai)GA(i ' ° ; °' ^ M .-+^ ( A7 ) 

where Ci\(L, M;m, fi) is the Clebsch-Gordan coefficient for adding the angular momentum quantum numbers (l,m) 
and (A,/x) and for forming (L,M). In our case, we have 

/, )m= ±i(fcz) = VM2i + i) y -^y2 e ^( k ' s )E(- i ) L ^( fcs ) 

x [Cn(L, 0; =f1, ±l)Gi(i, 0; 0, 0) - Cii(L, 0; ±1, Tl)C n (L, 0; 0, 0)] . (A8) 
Thus, the relevant coefficients are 



Cu(i + 1,0;0,0) = C,i(i + 1,0;±1,T1) = ±1 



2/+1' v ' ' ' ' V 2(2/ + 1)' 



\/^+T C ' l( ' - 1, 0; ±1, Ti) = W a^+i) • ( AQ ) 



Putting these together gives 



fi, m =±i(kz)^(- t ) l +^^±^- J -^L e -Tq ± (k,s)\j l+1 (ks)+j l - 1 (ks)] 

= H) i+ V^0 + 1)(2/ + 1) / ^e-g ± (k )S )^M. (A10) 



15 

Now, we get back to the observer's frame by applying the standard rotation operator, S(q), that takes the z-direction 
into fc. This gives 

/ lm (k) = D l m ^,(S(k))f lm ,{kz), (All) 

m'=±l 

where D , = (l,m'\S\l,m) is the matrix representation of the finite rotation of an initial state (l,m) into a final 
state (l,m ). We obtain 

a lm = -^ /(S E D^ASCm-^VWTTjWTT) I * e-iLCk,-)^. (A12) 

m' — ±1 

Finally, we calculate the angular power spectrum from (aj m a*/ m /) = CiSii>8 rn m' and obtain 

ft _ i£±i) ^e-/ ^-'"/^^(M), (A 13 , 

where P 9i is the power spectrum of q± defined by (2ir) 3 P q± (k)5 D (k — k') = (qj_(k)qj_(k')). Here, we have used the 
identity, 

/ Air 

d 2 k D l m:±1 (S(k))D l v l ±1 (S(k)) = -^—^S^Sa: (A14) 

The integral over k can be performed with Limber's approximation: when a function g(k, s) varies much more slowly 
than the spherical Bessel function, one can approximate the integral as 

^dkji(ka)ji(k^)g(k,a) « ^V^ ff (k = -,a). (A15) 



With this approximation, we finally obtain the desired formula for the kSZ power spectrum: 

Q = rarn^y f ds ^ ^(t = 1 1 'a, a) 
\ c J J s 2 et(s) 4 2 

This is Equation flBJ). 

CORRECTING FOR THE MISSING POWER IN SIMULATIONS 

The transverse momentum power spectrum at a given wavenumber, P qj _(k), receives contributions from the density 
and velocity auto/cross power spectra at various wavenumbers via Equation (|9]). As a result, P q± computed from a 
simulation with a fi nite box suffers fr om a loss of power due to the lack of modes whose wavelength is greater than 
the size of the box (jlliev et al J 1 20070 ). 

The missing power arises because we do not have P x n + g\ iX n+g\(k), P vv (k), or P x n+8),v(k) for k < fcbox = 27r/Zbox> 
where Zbox is the size of the box. In Equation ([9]) . this leads to the missing contributions in |k'| < fcbox and |k — 
k'| < fcbox- Estimating and correcting for the missing power thus requires the knowledge of the large-scale limit of 

Pxii+S),x(i+S) j Pwi an d P x ^ + s),v 

For the homogeneous reionization case, it is straightforward to recover the missing power, as the large-scale limits 
of P vv , P x (i + s),x(i+S) ( = X 2 Pss), and P x n + s),v{ = X.Psv) are precisely known by the cosmological linear perturbation 
theory. Using from the linear theory and the linear relation, P vv = (af /k) 2 P$$, we find that the missing-power- 
corrected momentum power spectrum from the A^-body simulation agrees precisely with the expected OV spectrum 
(see Figured]). Note that most of the missing power comes from Paa(|k — k'|)P„,;(A;') in fc' < fcb ox because of the 
relation, v(k) cx <5(k)/fc, in the large-scale limit. 

For the inhomogeneous reionization case, we do not have a precise way to calculate the ionized density power, 
P x (i+5) iX (i+5) , in the large-scale limit; however, we expect that the density field and the ionization field are reasonably 
flat at the scales larger than the box size, and correct for the missing bulk velocity of the box. Therefore, we expect 
that the term P x (i+gy x n+5\(\k — h'\)P vv (k') in fc' < fcb ox captures most of the missing power, as we have seen from 
the homogeneous reionization case above. With this approximation, the missing power in the inhomogencously ionized 
regime is given by 

f ri 3 k' o 

pM issi „ g(M) = __ (1 _^ )Px(i+5)x(i+5)(|k _ k/|)pTO(fc , ) _ (B1) 



In order to check the accuracy of Equation (|B1[) . wc compare the missing-power-corrected momentum power spectrum 
from the box of 114 ft, -1 Mpc (black solid line; denoted as L2) with that from a larger box of 425 Mpc (black 
dashed line; XL2) in Figure 21 We find a very good agreement between the two, confirming the robustness of our 
correction for the missing power. 
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